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We study several one dimensional step flow models. Numerical simulations show that the slope of 
the profile exhibits scaling in all cases. We apply a scaling ansatz to the various step flow models and 
investigate their long time evolution. This evolution is described in terms of a continuous step density 
function, which scales in time according to D{x,t) = F{xt~^^~'). The value of the scaling exponent 
7 depends on the mass transport mechanism. When steps exchange atoms with a global reservoir 
the value of 7 is 2. On the other hand, when the steps can only exchange atoms with neighboring 
terraces, 7 = 4. We compute the step density scaling function for three different profiles for both 
global and local exchange mechanisms. The computed density functions coincide with simulations 
of the discrete systems. These results are compared to those given by the continuum approach of 
Mullins. 
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I. INTRODUCTION 

The morphological evolution of crystalline surfaces has 
long been a major focus of attention in surface science. 
At low temperatures, below the roughening temperature 
of a low index facet plane, this is often dominated by the 
motion of surface steps, typically produced by miscuts or 
arising from dislocationsEl. Vicinal surfaces, created by a 
small miscut to the facet plane, have steps all of the same 
sign (either up or down steps). They offer a particularly 
simple testing ground for the study of step models of 
kinetic processes and their connections to physics both 
on atomic and macroscopic scales. 

In this paper, we study the kinetics of faceting and 
relaxation of vicinal surfaces using one-dimensional (ID) 
models of straight steps. We show that on long length 
and time scales the surface profiles exhibit a scaling be- 
havior in many cases. This general conclusion is not at 
all surprising; it agrees with the classic work by MullinsEl, 
who first investigated problems of this kind using a sim- 
ple continuum model. However, the details of the cal- 
culations we carried out and the nature of the scaling 
functions describing the surface profiles are rather dif- 
ferent. These differences arise because we take the con- 
tinuum limit of a physical model where steps play the 
fundamental role in the kinetic processes describing sur- 
face evolution. Even in the continuum limit there are 
features of the scaled profiles that reflect this underlying 
physics. There has been much discussion in the litera- 
ture about the .-derivation of continuum equations from 
discrete modelsu. Our work shows in detail how this can 
be done in some simple cases and illustrates some of the 
subtle issues that can arise. 

MuUinffl used a simple ID continuum model to de- 
scribe the growth of a single linear facet on a vicinal sur- 



face, assuming that the surface free energy of the complex 
surface away from the facet is an analytic function of the 
local surface slope. This assumption is valid in many 
cases since the vicinal surface itself is roughEl. However, 
since the existence of a facet is associated with the break- 
down of analyticity of the surface free energy as a func- 
tion of slope, this approach will fail whenever the local 
slope approaches zero, the slope of the flat facet. There 
is nothing in the usual continuum approach to prevent 
this from happening. Indeed, as noted below, in some 
cases the standard continuum model predicts that sur- 
face proflles will oscillate sufficiently about the slope of 
the vicinal surface to produce regions with negative local 
slope. This implies the creation of new anti-steps (steps 
of the opposite sign), an energetically costly process that 
does not occur in the step models we consider, even in the 
continuum limit. Furthermore, since at low temperatures 
steps play an essential role in determining surface kinetic 
processes, the effective kinetic coefficients in a continuum 
model must depend on the step density in ways that may 
seem hard to understand when viewing the system from 
a continuum perspective from the outset. 

To examine these issues in detail, we start with a ID 
model of discrete straight steps and describe surface mor- 
phological changes in terms of their motion. The equa- 
tions of motion for the individual steps reflect the mass 
transport mode and are derived in the standard waytJ, 
using a linear kinetics assumption based on the differ- 
ence in the chemical potentials for each step arising from 
step repulsions. By considering length scales large com- 
pared to the step spacing, we show that it is possible to 
take the continuum limit of these equations in a consis- 
tent way. The dynamic equation for the evolution of the 
local slope, or the step density function, is thus obtained 
systematically from the equations for individual step mo- 
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tion. The long time evolution of this dynamic equation is 
then investigated with a scaling ansatz. We applied this 
scheme to three different ID physical systems: recon- 
struction driven faceting, relaxation of an infinite bunch, 
and flattening of a groove. In all cases, the scaling func- 
tion is described by the same differential equation for 
the same mass transport mode; the solutions differ only 
because of different boundary conditions for the scaling 
functions. As a result, the values of the scaling exponents 
depend only on the mass transport mechanism, and are 
consistent with the predictions of MuUins' classical the- 
ory. However, the scaling functions themselves — the 
scaled slopes of the surface profile — differ from MuUins' 
results and are in excellent agreement with numerical so- 
lutions of the discrete equations. 

The paper is organized as follows. We introduce the ID 
step models in Sec. II. In Sec. Ill the step density func- 
tion is first defined. Then we introduce a scaling ansatz 
and derive a differential equation for the scaling func- 
tions. In Sec. IV, the properties of the scaling function 
are further investigated. In Sec. V, we check the valid- 
ity of the scaling analysis using three different step flow 
models. The shape of the scaling function is obtained 
by numerical integration of the differential equation. We 
also show that this scaling function coincides with the 
result of the simulations of the discrete step flow models. 
Our conclusions are given in Sec. VI. 



II. ONE DIMENSIONAL STEP FLOW MODELS 

Below the roughening temperature of a high symme- 
try orientation of the crystal, a vicinal surface consists of 
flat terraces separated by atomic steps. Ignoring islands 
and vacancies, the morphological evolution of the sur- 
face is a consequence of exchange of atoms between steps 
and their neighboring terraces resulting in motion of the 
steps. An analysis of surface evolution in terms of ste 
flow was first carried out by Burton, Cabrera and Fran^ 
and this has been generalized by many other authorsH. In 
what follows we describe the evolution of vicinal surfaces 
using a simple step fiow model. 

We consider two limiting channels for mass transport 
involving the terrace adatoms. In the first case the 
adatom mass flow on each terrace is local, and takes 
place by surface diffusion. In the second case the ter- 
race adatoms can easily exchange with a global reservoir, 
perhaps through direct hops to distant regions of the sur- 
face or by rapid exchange with the vapor. We refer to 
these limiting cases as the local (LEM) and global (GEM) 
exchange mechanisms, respectively. They correspond to 
the surface diffusion and evaporation-condensation mass 
transport mechanisms considered by MuUinaS. 



A. Local Mass Exchange Mechanism 

Consider an array of flat terraces separated by straight 
parallel steps with horizontal positions, The in- 

dex n grows in the direction of positive surface slope. 
These steps may absorb or emit atoms which then dif- 
fuse across the neighboring terraces. We ignore evapora- 
tion. Assuming attachment-detachment limited kinetics 
(i.e., that diffusion on terraces is very fast compared to 
the rate of attachment and detachment of atoms to and 
from step edges), the adatoms which diffuse on the nth 
terrace maintain a uniform chemical potential /i^ across 
the terracecl. 

We further assume that the flux of atoms at the two 
step edges bounding the nth terrace is determined by first 
order kinetics, characterized by a (for simplicity symmet- 
ric) attachment-detachment rate coefficient k: 



7^ 



= k 



(1) 



Here J+ and J,7 denote the flux from the lower and up- 
per neighboring terraces into the nth step respectively, 
/i* is the step chemical potential associated with ajddjng 
an adatom to the nth step. In the case of elasticQ'El or 
entropic repulsive interactions between steps, it is well 
knowuthat the step chemical potential then takes the 
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(Xn ^n — 1)^ ^n)"^ 



(2) 



where g is the strength of the repulsive interactions. 

Next, we assume that diffusion processes are fast com- 
pared with the motion of steps. Within this quasi-static 
approximation, the density of adatoms on the terraces 
reaches a steady state for each step configuration. In 
this steady state we have + J,7+i = for any n, and 
therefore the free energy associated with the addition of 
an adatom on the nth terrace, /i^, takes the form 



(3) 



Combining mass conservation at the nth step with Eqs. 
(|l|) and @), we obtain the following expression for the 
step velocity: 



c,i a ( J,^ + J,j 



(4) 



with a denoting the lattice constant of the crystal. 



B. Global Mass Exchange Mechanism 



As with the LEM, we assume linear kinetics; i.e., the 
flux of atoms from the reservoir to the nth step is pro- 
portional to the difference between chemical potentials of 
the reservoir and the step: 
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(5) 



Here fires is the reservoir chemical potential, which we 
choose as the chemical potential of a flat surface; i.e., 
I^res = 0. The velocity of the nth step in this case is 
simply 



-a Jn — a kfi:^ 



(6) 



Although the exchange rates in the LEM and GEM cases 
are different we use the symbol k to denote both of them. 
It should be clear from the context which exchange rate 
we refer to. 



III. SCALING ANALYSIS AND CONTINUUM 
MODELS 



that the scaling ansatz Eq. (|^) holds, we now change 
variables to 9 = t^l^ and ^„ = a;„0~^, and transform Eq. 
into an equation for the scaling function F evaluated 

atC= (6^+en+i)/2: 



dF 



+ F'^Q-^ (i„+i - i„) = , 
(10) 

In going from Eq. (||) to Eq. (|lo|) we have used the fact 
that dD jdt = —D'^{xn+i — in) and that dx/dt = {in + 
in+i)/2. Xm Xn+i thcmsclves are now expressed in term 
of the ^„ 's using Eq. (^ or (||) depending on the exchange 
mechanism. 

Let us also rewrite Eq. (^) in terms of 6*, F, and the 
6;s: 



Simulations of the LEM and GEM step flow modelsB 
suggest that the behavior of these systems at long length 
and time scales can be described in terms of a step den- 
sity function D (i.e., the inverse step separation) that is 
a continuous function of both position and time. Gen- 
erally speaking, continuum descriptions of step systems 
should be valid when every typical surface feature con- 
sists of very many steps. In the cases we study here, the 
simulations indicate that the step density function scales 
in time according to 



D{x,t) = F{xt-^'^) 



(7) 



with a positive exponent 7. Thus, the number of steps in 
every surface feature grows with time as t^/'^ . We should 
therefore be able to accurately describe the evolution of 
the system in terms of a continuum model in the long 
time limit. 

Eq. (|^) makes the even stronger assertion that the sur- 
face features have a self-similar shape during their evo- 
lution, and can be related by a proper rescaling of time 
and distance. This scenario is similar to the scaling ex- 
hibited by a, dec aying crystalline cone, which two of us 
have studiedli30. In this section we carry out a scaling 
analysis, similar to the one in Rcfs. |lO| and |l^, to obtain 
the scaling exponents 7 and the differential equation for 
the scaling function F. We also study the effects of the 
different mass exchange mechanisms. 

We start by defining the step density function in the 
middle of the terraces: 



D ( ^"^^"+\ t) = 



(8) 



Assuming continuity, the full time derivative of the step 
density is given by: 



dD 
~dt 



dD 

'at 



dD dx 
dx dt 



(9) 



Eq. (|^) can be evaluated in the middle of the terrace 
between two steps (i.e., at a; = (a;„ + a;„+i)/2). Assuming 



n-l-l 



en 



F[(^„+i+Cn)/2] 



(11) 



According to this, the difference between successive ^n's 
is of order wherever F does not vanish. In the large 
9 (long time) limit these differences become vanishingly 
small. The differences ^n+k — £, will also be small as long 
as k is finite. We can therefore take the continuum limit 
in the variable e, consistent with our original supposition. 

To this end, we evaluate the function F at the position 
(6i+fc + eri+fe+i)/2 by using its Taylor expansion 



S,n+k + Cn+fe+l 



9- 



1 d^f (0 ^ in+k+in+k+1 , 



m=0 



TO 



de 



Next, we expand 



>kmi 



(12) 



(13) 



and insert this into Eq. (|12|). By equating terms of the 
same order in 9^^ on both sides of Eq. (|l2|), we can find 
the coefficients (j)kin for any desired values of k and to. 
These coefficients involve the function F and its deriva- 
tives evaluated at ^ = (^„ -I- ^„_|_i)/2. 

Having found the expansion coefficients (pkm we now 
return to Eq. (|l^) . This equation depends on the veloci- 
ties Xn and Xn+i which in turn depend on ^71-2 ■•■ ^"+3 
or ^n-i ••■ £n+2 in the LEM or GEM cases respectively. 
Using Eq. ( [l3| ) we expand Eq. (^ in the LEM case with 
the following result: 



9''-^a'^kg— 



1 d2 3_p2 



e dF 
7 d^ 



r-^) = . 

(14) 



Consider Eq. (|14D. It involves different powers of the 
scaled time 9 and cannot be satisfied at all times (unless 
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F is trivially independent of ^). However in the 9 
oo (long time) limit Eq. can be satisfied exactly. 
This can be achieved by setting the value of the scaling 
exponent to be 7 = 4 and then requiring the first and 
second terms to cancel each other. Neglecting the small 
third term we are left with the long time limit of the 
LEM differential equation 



Aa kg- 



1 3F^ 



\FdC 2 



.dF 
di 







(15) 



This equation determines the scaling function F and is 
exact only in the long time limit. In this limit it is equiva- 
lent to the attachment/detachment limited scaling equa- 
tion suggested by Liu jSt all3. It also agrees with Nozieres' 
continuum treatmentllj of attachment/detachment lim- 
ited kinetics provided one assumes scaling. Thus, in the 
long time limit, our analysis confirms the validity of pre- 
vious continuum treatments. However at any finite time 
there are corrections of order to the scaling solution 
Dix,t) = F{i). 

A similar procedure can be applied in the GEM case. 
It leads to the scaling exponent 7 = 2 and the following 
equation for the scaling function: 



, 2, d^F^ JF 







(16) 



As in the LEM case, the leading correction to the density 
fmiction decays as 9^^. Like MuUins, we obtain a fourth 
order equation with LEM and a second order equation 
with GEM. However, these equations have a more com- 
plicated form than those arising from Mullins' continuum 
model and will have different solutions. 

To conclude this part let us emphasize a few points 
regarding the above analysis. Eqs. ( p5| ) and (|l^) to- 
gether with the scaling ansatz Eq. (^ and the values of 
7 constitute a continuum model for the surface dynam- 
ics. This model was derived directly from the discrete 
step system and is exact in the long time limit. In addi- 
tion, the above scaling analysis is robust in the following 
sense. The values of the scaling exponent 7 are not sen- 
sitive to the exact nature of the step-step interactions in 
the discrete model. In Appendix A we show that for a 
general interaction 7 = 4 and 7 = 2 in the LEM and 
GEM cases respectively. However, the differential equa- 
tions for the scaling functions do depend on the form of 
the interaction. 



IV. PROPERTIES OF THE SCALING FUNCTION 

Here we study some properties of our continuum 
model. Our purpose is to derive several relations for the 
scaling function, which hold in general. These relations 
are useful in the derivation of the boundary conditions 
necessary in order to solve Eqs. ( p^ and ( p^ ) for the 
scaling functions of various systems. 



First, note that the scaling function must have a finite 
limit, Foo, at infinity. Otherwise the step density there 
would change infinitely fast. We choose the unit of length 
so that Foo — 1. We also choose the unit of time so that 
ia^kg = 1 in Eqs. (|l|) and (|l|). 

Next, we investigate the time dependence of the vol- 
ume and the total number of steps in the system. It turns 
out that these quantities can be calculated directly from 
the differential equations. The change in the volume of 
the system in the positive x half of space during the time 
interval from zero to t is given by 



AV 



[h{x, t) — h(x, 0)] dx 



(17) 



h{x,t) ~ D{x,t)dx is the profile height measured in 
units of the lattice constant. 

Integrating by parts and changing to scaling variables 
we obtain the equation 



AV = r 



iFoo-F)^dC 



(18) 



where we have ignored the surface term assuming no evo- 
lution occurs infinitely far from the origin. 

Let us calculate the last integral in the LEM case. Ac- 
cording to Eq. (nSl) 



^de \Fde 2 ' ^ 



ill rip 

d^ 



(19) 



where and ^f, are boundaries of a region where Eq. 
( [Tsl ) is valid. Integrating each term by parts we find that 



AV 



LEM 



e {Foo -F)+^ 



1 d^ 3F^ 



d^ yFd^"^ 2 

{F^ - F) e^e ■ 



1 d^ 3F^ 



Fde 



(20) 



Similarly, we can calculate the volume change in the 
GEM case using Eq. dl^) and obtain the equation 



AV( 



GEM — 



dF^ 

e {Fo.-F)-e-^+F^ 

d^ 



+9^ {F^-F)^d^ 



(21) 



A similar treatment can be applied to calculate changes 
in the number of steps in the positive x half of the system, 
and we find 



/•oo 

AN = 9 {F-F^)d^ 
Jo 



(22) 



Again we can evaluate this integral by using Eqs. ( [T5| ) 
and ( n™. The results are 
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A7V, 



LEM 



1 (P 



in the LEM case and 



{F-Foo)di 



(23) 



GEM 



dF^' 

HF-F^) + ^ 



+e [ {F-Fo,)dc 



(24) 



in the GEM case. 

Finally we study the behavior of the scaling functions 
near regions of zero step density. We recall that our scal- 
ing analysis is valid only in regions of space where the step 
density does not vanish (see Eq. (0)). Points of vanish- 
ing step density should therefore be treated separately. 
Assume that is such a point for which F{£^o) — 0. Ex- 
panding the scaling function in powers of ^ — in the 
vicinity of and using Eqs. (|l^) and (|l6|), we find that 
in the LEM case 



i/2 



and in the GEM case 

C30 



(25) 



(26) 



n=l 



Thus a point of vanishing step density is a singular point 
of the scaling function at which all its derivatives diverge. 



V. EXAMPLES 

To check the validity of our scaling analysis, we con- 
sider three different step flow models, which according to 
numerical simulations obey scaling under both exchange 
mechanisms. All three systems consist of straight par- 
allel and initially equidistant steps, but differ in their 
boundary conditions. 



configuration where all the terraces except one have the 
same width w < Wc- (For simplicity, we do not permit 
reconstruction on other terraces, even if their widths ex- 
ceeds Wc during the surface evolution. The possibility of 
such "induced nucleation" on other terraces is discussed 
in Refs. ^ and |l4[) The zeroth terrace (between xq and 
Xi in Fig. |l| (a)) is of different width, larger than Wc 
and is reconstructed. This reconstructed terrace tends 
to become even wider, and this is reflected as a shift of 
magnitude e in the chemical potentials of the two steps: 



Mo = 2.g 



1 



1 



{xo - x-i)^ {xi - xq)^ 
1 1 



(Xi-Xo)^ ix2-Xiy 



(27) 



As a result, the step at xq propagates to the left, while 
the one at xi propagates to the right. In the long time 
limit the step density of the system obeys scaling. In 
particular, the size of the facet at the origin grows as 
t^^f in both LEM and GEM cases, with 7 = 4 in the 
LEM case and 7 = 2 in the GEM case. 




FIG. 1. The three systems considered in this section: (a) 
reconstruction driven faceting, (b) relaxation of an infinite 
bunch and (c) flattening of a groove. 



A. Reconstruction Driven Faceting 

The first example is that of reconstruction driven 
faceting studied by Jeong and Weeks in Refs. || and Q It 
is mathematically equivalent to the modeljof facet growth 
during thermal etching studied by MuUinsH. Within their 
model, surface reconstruction that lowers the free energy 
can nucleate only on terraces of width larger than a crit- 
ical width, Wc- Therefore, terraces of width w > Wc have 
a lower free energy than those of width w < Wc- We con- 
sider the evolution of such a step system starting from a 



In order to compare our scaling analysis to simulation 
results we have to solve Eqs. (|l^) and ( p^ ) with the rele- 
vant boundary conditions. Since the system is symmetric 
about the origin it is sufficient to solve the scaling func- 
tion F for positive ^. In addition, our expansion in the 
small parameter is valid only in regions where F does 
not vanish (see Eq. (pd|)). This requirement is violated 
on the diverging facet around the origin. Therefore, F 
obeys Eqs. (^ or (^ only for ^ > ^1, where ^1 is the 
scaled position of the first step, and F(^) = for ^ < Ci- 
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1. Local Mass Exchange Mechanism 



.1/3 



(32) 



The first boundary condition is set by our choice of the 
value of F at infinity, namely f oo = 1 • 

Next, consider the number of steps, which is a con- 
served quantity in our step model. This is an important 
difference from the continuum approach of MuUins. Us- 
ing Eq. ( p3| ) with £,a = £.b = oo and the condition 
F(^) = for ^ < ^1, we find that in the LEM case 



F'"(a) = ^ - 



2F'F" 
F2" ~ 



3 



(28) 



where primes denote derivatives with respect to ^. 

Two additional boundary conditions can be found by 
investigating the velocity of the first step. According 
to our scaling ansatz the position of the first step is 
a^i = ^^1- In the long time scaling limit the velocity 
of this step goes to zero as 9~^. Using Eq. (||) together 
with the symmetry of the system we obtain the following 
expression for ii: 



1 



1 



(2xi)3 (X2-Xi)3 (a;3-X2)3 



3e 



(29) 



Rewriting Eq. (|29| ) in terms of the scaling variables 
and expanding in we find that 



Xl 



F^) 



iFF'e- 



o{e- 



(30) 



where F is evaluated at (^i + £,2)/2. In order for xi to 
vanish as we must have F[{^i + C2)/2] — f^l'^ and 
F'\{i\ + 6)72] = 0. Terms of order Q-"^ on the r.h.s. also 
have to vanish, but they include corrections to scaling, 
which we ignore in this work. In the long time limit the 
difference between and (^i -l-^2)/2 is negligible and we 
arrive at the boundary conditions 



nil) 



and F'(^i) = 



(31) 



We solved Eq. (y_5|) numerically by applying the three 
boundary conditions at ^ = (Eqs. ( p8| ) and (|3l|)), and 
by tuning -F"(^i) and itself to satisfy the boundary 
condition at infinity. In the upper half of Fig. |^ we com- 
pare the resulting solution with scaled density functions 
taken from numerical simulations of the discrete model. 
The agreement is quite impressive. These results should 
be compared with Fig. 8 in Ref. |^, which predicts anti- 
step formation for large values of the slope parameters 
m/n. 



2. Global Mass Exchange Mechamsm 

Like in the LEM case, the requirement that the velocity 
of the first step decays in time results in the boundary 
condition 



Again is the scaled position of the first step. 

Another boundary condition is derived from the con- 
servation of the number of steps. Eq. ( p4| ) with = ^i, 
^6 = cx) and the condition F(^) = for ^ < implies 
that 



F'iii) = - 



(33) 



These two conditions are sufficient for solving Eq. ( p^ ) 
when the value of is known. 

Finally we tune to satisfy the boundary condition at 
infinity. The resulting solution and its comparison with 
simulation data are shown in the lower half of Fig. |^. 
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FIG. 2. Solution of the reconstruction driven faceting scal- 
ing function in the LEM and GEM cases (solid lines) com- 
pared with scaled density functions from numerical simula- 
tions with 6 = 8. Different symbols represent density func- 
tions at different times. 



B. Relaxation of an Infinite Bunch 

In the second example we study the relaxation of an 
infinite bunch of steps. The initial step configuration 
(Fig. |l| (b)) consists of an infinite facet at a; < 0, in 
contact with an infinite array of uniformly spaced steps 
(at X > 0). The first step has a single neighbor, and 
therefore its chemical potential is 



Ml 



-2ff 



{X2 - XlY 



(34) 



In the LEM case there is a complication, since the first 
step may exchange adatoms with the infinite facet on its 
left. Hence we have to specify /j,q, the adatom chemical 
potential on the facet. We assume here that /Xg = /xf, 
neglecting any exchange of atoms between the first step 
and the facet. As a result the volume is conserved in the 
LEM case. 



6 



Simulations of this system in both the LEM and GEM 
cases suggest that the system exhibit scahng with the 
origin of the scaled position at the initial position of the 
first step. The leftmost steps from the bunch move to 
the left into the facet due to the repulsive interactions. 
The first step recedes in time and its position scales as 
xi — i^/^. At the same time the separation between 
the first steps grows. For ^ < ^i, the scaled position of 
the first step, the step density always vanishes, and we 
have to find the scaling function only for < ^ < cxd. 



1. Local Mass Exchange Mechanism 

The velocity of the first step vanishes in the scaling 
limit. We can therefore evaluate F{^i) by expanding xi 
in 6~^, and requiring that the zeroth order term vanish. 
We find that the zeroth order term in Xi is proportional 
to F^{£_i), which implies that i^(Ci) = 0. 

Two additional boundary conditions can be derived 
from the conservation of volume and the number of steps. 
Using Eqs. (^0|) and (^3|) to calculate changes in the vol- 
ume and number of steps in the positive part of the ^ axis 
together with equivalent equations for ^ < 0, we find that 



di \Fde 2 



1 3F2 



Fd£_^ 



and 



d_ 



]_d^3F^ 
Fd^~ 



= 



, (35) 



(36) 



As mentioned above, a zero of the scaling function is a 
singular point at which all derivatives of F diverge. Since 
F{^i) = 0, the limit ^ ^ Ci in the last two boundary 
conditions should be taken with care. This can be done 
by considering the power series (^5|). The differential 
equation (|l5|) and the boundary conditions -F'(Ci) = 0, 
( Pq ) and ( |36^ impose connections between the coefficients 
of the expansion and leave only bi as a free parameter. 

We now use the following procedure to compute F. For 
given values of bi and we approximate F at £^i+S^ us- 
ing the series expansion ( p5| ) with suitable truncation. At 
^1 -|-(5^ the derivatives of F are finite, and we solve F from 
there by numerical integration. We adjust the values of 
bi and in order to satisfy the boundary condition at 
infinity, thus obtaining the scaling function -F(^) in the 
LEM case. In the upper half of Fig. || we compare this 
solution with simulation data. 



2. Global Mass Exchange Mechanism 

Since xi oc — (a;2 — xi)~^, the requirement that the 
velocity of the first step vanishes in the scaling limit leads 
to the boundary condition 



F{^i) = . 

To derive another boundary condition we use Eq. 
and impose conservation of the number of steps, 
following relation is thus obtained: 



dF^ 







(37) 
The 



(38) 



This implies that all the coefficients in the series expan- 
sion (|2^) diverge except for bi which vanishes. Thus the 
convergence radius of ( p6|) is zero and it is difhcult to 
calculate the scaling function with the numerical proce- 
dure used in the LEM case. Nevertheless, taking a small 
enough value of bi and tuning to satisfy the bound- 
ary condition at infinity, we were able to calculate F 
approximately. In the lower half of Fig. ^ we compare 
this approximation with simulation data. 
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FIG. 3. Solution of the infinite bunch scaling function in 
the LEM and GEM cases (solid lines) compared with scaled 
density functions from numerical simulations of the discrete 
model. Different symbols represent density functions at dif- 
ferent times. 



C. Flattening Of A Groove 

Our last example is the flattening of a groove cut in 
the crystal surface. The initial conflguration (Fig. Q (c)) 
consists of two infinite step bunches with steps of oppo- 
site signs that meet at the origin. Step repulsion within 
each bunch pushes the bottom step and anti-step towards 
each other until they collide and annihilate. We assume 
here that steps of opposite sign do not interact. Thus 
the chemical potential of the two bottom steps includes 
interaction with only one neighboring step, namely 

fil = , .3 . (39) 

{X2 - XiY 

When the first steps annihilate we relabel the remaining 
steps so that the positions of the bottom step and anti- 
step are always Xi and xq. Our system is symmetric with 



7 



respect to the origin, and it is sufficient to consider only 
the positive part of the x axis. Symmetry also excludes 
any flux of adatoms between the two bunches, and in 
the LEM case this implies that the volumes of the two 
bunches are conserved separately. 



1. Local Mass Exchange Mechanism 

The groove flattening problem is different from the two 
previous examples, where there was a well defined facet 
with its edge at the position of the first step. At the 
facet edge Eq. (|l^) becomes invalid due to the vanishing 
step density. In the groove example, steps annihilate at 
the origin, and there are steps in all regions of space. 
However, there could still be points where the density of 
steps vanishes if the position of the second step diverges 
in the scaling limit. Let us denote the scaled position of 
the facet edge by i^*. Eq. is vahd only for ^ > £,*, 
and the value of ^* is unknown a priori. There are two 
possible, qualitatively different situations: 1) ^* > 0, i.e., 
there is a plateau at the bottom of the groove which 
grow as t^^^. 2) ^* ~ 0. There could still be a diverging 
plateau, but it must grow more slowly than t^^'^ . In what 
follows we rule out the first possibility and show that 

r = 0. 

According to Eq. (p3|), the number of step annihilation 
events grows with time as t^/''. Denoting by tn the time 
of the nth annihilation event, we see that tn ~ n'^ and 
tn+i — tn ^ n^. We now show that if £,* > 0, the time 
interval between annihilation events is larger than n^. 
Just before the nth annihilation event, the velocity of 
the first step (which must be negative) is given by 



lim xi 



1 



(X3 - X2 y x-^ 



If > 0, X2 is of order tl/"^ ~ n. Eq. (^) imphes that 
X3 ~ X2 > 2^^/^X2, and therefore the distance 2:3 — X2 is 
also of order n. 

After the nth annihilation event we relabel the steps. 
X2 becomes xi , X3 becomes X2 and so on. Now the dis- 
tance between the first two steps, X2 — xi is of order n. 
The velocity of the new first step, which is maximal (in 
absolute value) at this time, obeys the following inequal- 
ity: 



(40) 



lim < 



1 



2{x2-xiy 



0{r 



(41) 



The step must cross a distance of order n until it anni- 
hilates and therefore tn+i — tn is at least of order n*, in 



tr. 



derived 



contradiction with the relation 
above. Hence ^* = 0. 

In order to solve Eq. (p^), we now find two boundary 
conditions at ^ = 0. Consider first the quantity X3 — X2- 
If X2 diverges with time, x^ — X2 must also diverge in 
order for xi to be negative (see Eq. (40)). If X2 does 



not diverge in the scaling limit, x^ — X2 must still diverge 
in order for the time interval between two consecutive 
annihilation events to diverge in the scaling limit. Thus, 
X3 — X2 diverges and the scaling function vanishes at the 
origin: 



F(0) = 



(42) 



Another boundary condition is derived from volume 
conservation in the x > half of space. From Eq. ( [20| ) 
we obtain the condition 



1 3f ^ 



1 iF^ 
Fd^~ 



= 



(43) 



Requiring expansion (g5|) to satisfy Eqs. (|15|) and (|42 
we are left with two free expansion coefficients, which 
we tune in order to satisfy the boundary condition at 
infinity. The resulting solution compared to simulation 
data is shown in the upper half of Fig. |. 




FIG. 4. Solution of the groove flattening scaling function in 
the LEM and GEM cases (solid lines) compared with scaled 
density functions from numerical simulations of the discrete 
models. Different symbols represent density functions at dif- 
ferent times. 



2. Global Mass Exchange Mechamsm 

As in the LEM case, we have to find the value of 
From considerations similar to those used in the LEM 
case it can be shown that steps annihilate fast enough 
for scaling to occur only if ^* =0. 

In the GEM case the volume is not conserved. How- 
ever, by summing up the velocities of all the steps we can 
calculate the rate of change of the volume: 



dy_ 

'dt 



n=l 



— lim 

2 Af— 00 



{XN - XN-lY 



(44) 



The r.h.s. of Eq. 

calculated from Eq. 



is F^/2, while the l.h.s. can be 
y) with = 0, = 00 and 9^ = 



8 



t. This calculation combined with Eq. 
boundary condition 



3 

oo 7 



which implies 



FiO) = 0. 



|) leads to the 
(45) 

(46) 



To calculate the scaling function we use the series ex- 
pansion ( p^ ) to approximate F near the singular point 
and then integrate Eq. ( p^ numerically from a point 
where F is finite. We tune the expansion coefhcient bi 
to satisfy the boundary condition at infinity. The result- 
ing solution compared to simulation data is shown in the 
lower half of Fig. |^. 



VI. SUMMARY AND DISCUSSION 

We studied ID step-flow models for the kinetics of 
faceting and relaxation and showed that the surface pro- 
file exhibits scaling behavior. The value of the scaling 
exponents was determined by the mass transport mech- 
anism. The scaling functions in all cases considered here 
were described by the same differential equation for the 
same mass transport mode. These scaling functions dif- 
fer from those predicted by MuUins' continuum theory, 
which does not explicitly consider the existence of steps. 
However, the scaling exponents 7 agree with MuUins' 
classical theory. This can be understood as follows. The 
systems considered here are well below the roughening 
temperature of the (low-index) singular surface but large 
parts of the system are rough, with nonzero macroscopic 
slope and a differentiable free energy as a function of ori- 
entation. Therefore, if the system shows scaling behavior 
at all, the whole system (including the singular region) 
must evolve with the same time dependence as the non- 
singular part, which is accurately described by MuUins' 
model. 

Here we only considered attachment/detachment lim- 
ited kinetics for local mass exchange. We have carried out 
a similar scaling analysis for diffusion limited kinetics. 
Diffusion limited kinetics is also a local mass exchange 
mode and hence the scaling exponents are the same as for 
attachment/detachment limited kinetics. However, the 
differential equation for the diffusion limited scaling func- 
tion is different from the attachment/detachment limited 
case. For diffusion limited kinetics we found equations, 
which under the scaling assumption, are equivalent to the 
continuum models proposed in Refs. 0,^. Our results are 
valid even for finite step permeability at sufficiently long 
times. 

Another issue that we would like to investigate further 
is the condition for scaling behavior. When a step profile 
shows a scaling behavior, the scaling exponent should 
depend on the mass transport mechanism but not the 
driving force of the evolution at the boundary. However, 



very strong driving forces could cause a piling up of steps 
and destroy the scaling behavior. What kind of driving 
force gives rise to a scaling behavior? For example, for 
the flattening of a groove, we used a contact interaction 
between steps with opposite signs. What kind of inter- 
action between the steps of opposite signs in the middle, 
in general, results in a profile that obeys scaling? Does 
it depend on the mass transport mechanism? We have 
examined some artificial examples (not described here) 
of interactions between the steps of opposite signs that 
can destroy the scaling behavior but do not know yet the 
general form of the driving forces that admit a scaling 
ansatz. 

One can also consider more general interactions be- 
tween neighboring steps in the non-singular region. Al- 
though our scaling analysis is valid for a general step- 
step interaction (see Appendix A), we only carried out 
simulations to verify that scaling indeed occurs in the 
case of simple entropic and elastic step-step interactions. 
What are the interactions that support scaling? What 
is the general relationship between the boundary condi- 
tions and the step interactions that is consistent with a 
scaling ansatz? If these questions can be answered, we 
may be able to tell in advance which surface systems will 
show dynamic scaling behavior. 
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APPENDIX A: 

In this appendix we give the results of a scaling analysis 
with a general step chemical potential formula. We start 
by replacing Eq. (^ with 



U {Xn - Xn-l) - U {Xr. 



(Al) 



with U{x) a general analytic function of x. Such a 
formula is consistent with interactions between nearest 
neighbor steps. 

Next we use this formula together with the step ve- 
locities (|) and d) to rewrite Eq. (|l^ in the LEM and 
GEM cases. Changing the step velocities does not alter 
the expansion ([l^) since this expression is general and 
depends only on the scaling function F. We can thus use 
expansion (|l^) to expand the new Eq. (|lO|) in powers of 
0-^. We find that in the LEM case 



2 de 

7 




(A2) 



This implies that in the LEM case 7 = 4 and the differ- 
ential equation for the scaling function is 
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2a^k 



^_d_ ( 1 dU{^) ' 
Fd^ [f d^ 



.dF „ 



(A3) 



In the GEM we find that 
which implies 7 = 2 and 



(A4) 



2a^k- ^-^^ 



(A5) 
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